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A SPHERE  USING  AN  IMMERSED  BOUNDARY  METHOD 
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School  of  Mechanical  and  Aerospace  Engineering 
Seoul  National  University,  Seoul  151-742,  Korea 


Abstract.  Large  eddy  simulation  of  turbulent  flow  over  a sphere  is  con- 
ducted at  two  different  Reynolds  numbers  (Re  = 3700  and  101 * * 4)  based  on 
the  free-stream  velocity  and  the  sphere  diameter,  in  order  to  investigate 
the  flow  characteristics  over  the  sphere.  The  present  numerical  method  is 
based  on  a newly  developed  immersed  boundary  method  in  a cylindrical  co- 
ordinate (Kim  et  al.  2001).  Also,  a hybrid  numerical  method  of  eliminating 
2-A  waves  in  laminar  accelerating  flow  region  and  also  avoiding  numerical 
dissipation  in  turbulent  flow  region  is  newly  proposed  for  large  eddy  simu- 
lation of  turbulent  flow  over  a bluff  body.  At  Re  = 3700,  the  shear  layer  is 
elongated  in  the  streamwise  direction  to  form  a cylindrical  vortex  sheet  and 
its  instability  begins  to  appear  at  x « 2d.  Also,  the  flow  behind  the  sphere 
is  nearly  laminar  at  x < d and  contains  few  vortices.  On  the  other  hand, 
at  Re  = 104,  the  shear  layer  instability  occurs  right  behind  the  sphere  in 
the  form  of  vortex  tube  and  the  flow  becomes  turbulent  in  the  near  wake. 
Therefore,  at  Re  = 10,  the  size  of  the  recirculation  region  is  smaller  and 
the  wake  recovers  more  quickly  than  at  Re  = 3700. 


1.  Introduction 

Flow  behind  a sphere  shows  completely  unsteady  three-dimensional  phe- 

nomena even  at  low  Reynolds  numbers  (Sakamoto  & Haniu  1990;  Johnson 

& Patel  1999;  Mittal  1999).  Therefore,  flow  over  a sphere  has  not  been  well 
understood  as  compared  to  that  over  a circular  cylinder,  even  though  the 
sphere  geometry  is  simple.  Previous  studies  on  flow  over  a sphere  have  been 
mostly  based  on  experimental  methods,  where  vortical  structures  behind 
the  sphere  were  presented  using  flow  visualization  and  the  high  and  low  fre- 
quency characteristics  in  the  shear  layer,  respectively,  associated  with  the 
shear  layer  instability  and  the  wake  instability  were  investigated  (Sakamoto 
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& Haniu  1990).  When  the  Reynolds  number  exceeds  about  800,  the  wake 
flow  behind  a sphere  changes  from  laminar  to  turbulent  flow.  In  recent 
years,  a few  numerical  studies  have  been  made  on  turbulent  flow  over  a 
sphere  (Tomboulides  1993;  Constantinescu  & Squires  2000).  However,  the 
detailed  quantitative  information  on  the  turbulence  statistics  behind  the 
sphere  is  still  limited  in  the  literature. 

The  objective  of  the  present  study  is  to  investigate  the  characteristics 
of  turbulent  flow  over  a sphere  using  large  eddy  simulation  at  two  different 
Reynolds  numbers  of  3700  and  104  based  on  the  free-stream  velocity 
and  the  sphere  diameter  d. 

2.  Numerical  Method 

The  present  numerical  method  is  based  on  a newly  developed  immersed 
boundary  method  in  a cylindrical  coordinate  (Kim  et  al.  2001),  where  the 
sphere  in  the  flow  is  treated  as  momentum  forcing  in  the  Navier-Stokes 
equations.  Therefore,  flow  over  the  sphere  can  be  easily  handled  with  or- 
thogonal (cylindrical)  grids  that  do  not  coincide  with  the  body  surface. 
Numerical  simulations  are  performed  at  two  different  Reynolds  numbers 
of  Re  = Uood/v  = 3700  and  104.  The  computational  domain  used  is 
— 15d  < x < 15 d,  0 < r < 15 d,  0 < 9 < 2n,  and  the  number  of  grid 
points  is  577(x)  x 141  (r)  x 40(0).  Here,  x,  r and  0,  respectively,  denote  the 
streamwise,  radial  and  azimuthal  directions.  A Cartesian  coordinate  system 
(a;,  y,  z)  is  also  defined  in  order  to  present  the  drag  and  lift  forces,  where 
the  y and  z coordinates  coincide  with  0 = 0°  and  0 = 90°,  respectively. 

One  of  the  important  issues  for  large  eddy  simulation  of  turbulent  flow 
over  a bluff  body  at  high  Reynolds  number  is  that  the  second-order  cen- 
tral difference  scheme  produces  serious  dispersion  errors  for  convection- 
dominated  flow,  and  thus  2-A  waves  appear  in  the  region  of  velocity  accel- 
eration (for  example,  at  front  surface  of  a sphere).  To  resolve  this  problem, 
two  different  approaches  have  been  taken:  One  is  to  apply  a filtering  pro- 
cedure which  corresponds  to  adding  numerical  dissipation  in  an  ad  hoc 
manner,  and  the  other  is  to  use  upwind  or  upwind-biased  schemes.  In  these 
approaches,  however,  the  numerical  dissipation  becomes  larger  than  the 
dissipation  from  the  subgrid-scale  viscosity  (Beaudan  & Moin  1994). 

In  this  study,  we  propose  a hybrid  approach  for  the  spatial  discretiza- 
tion of  the  convection  terms  for  flow  over  a bluff  body,  where  a high- 
order  upwind  scheme  is  used  in  the  region  of  laminar  accelerating  flow  and 
the  second-order  central  difference  scheme  is  used  elsewhere.  The  diffusion 
terms  are  discretized  with  the  second-order  central  difference  scheme.  With 
this  approach,  we  can  avoid  possible  contamination  of  solutions  resulting 
from  the  application  of  the  dissipative  schemes  to  the  whole  flow  field. 
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(a)  (b) 

Figure  1.  Contours  of  the  instantaneous  radial  velocity  in  front  of  the  sphere  at  Re  = 
3700:  (a)  second-order  central  difference  scheme  (CD);  (b)  hybrid  scheme  (CUDZ3+CD). 


For  the  high-order  upwind  scheme,  we  use  the  third-order  compact  upwind 
scheme  (CUDZ3)  by  Zhong  (1998),  which  is  stable  and  less  dissipative  than 
the  straightforward  compact  upwind  scheme  using  an  upwind-biased  grid 
stencil  (Tolstykh  1994).  For  the  flow  over  a sphere,  we  apply  CUDZ3  to  the 
laminar  accelerating  flow  region  ( x < 0 and  r < d). 

Figure  1 shows  the  contours  of  the  instantaneous  radial  velocity  in  front 
of  the  sphere  at  Re  = 3700  using  two  different  numerical  schemes.  Figure 
1 (a)  is  obtained  from  the  second-order  central  difference  scheme  (CD) 
applied  to  the  whole  flow  field,  whereas  figure  1 (b)  is  from  the  hybrid 
scheme  proposed  in  the  present  study.  It  is  clear  that  the  present  hybrid 
scheme  effectively  removes  the  wiggles  appeared  in  figure  1 (a).  Therefore, 
all  the  following  results  presented  in  this  paper  are  obtained  using  the 
hybrid  scheme. 

3.  Results 

The  simulation  results  are  summarized  in  table  1,  where  time-averaged 
drag  coefficient  ( Cd ),  base  pressure  coefficient  ( CPb ),  separation  angle  (ds), 
recirculation  bubble  size  (L),  and  Strouhal  number  (St)  associated  with 
the  wake  instability  are  presented  together  with  experimental  and  compu- 
tational data  available.  Here,  LES  and  DES  denote  large  eddy  simulation 
and  detached  eddy  simulation,  respectively,  and  DES  is  a hybrid  approach 
reducing  to  Reynolds-averaged  Navier-Stokes  (RANS)  near  the  wall  and 
LES  away  from  the  wall.  The  recirculation  bubble  size  (L)  is  the  distance 
from  the  base  of  the  sphere  to  the  point  where  the  time-averaged  stream- 
wise  velocity  is  zero.  The  drag  coefficient,  base  pressure  coefficient,  Strouhal 
number  and  separation  angle  are  in  reasonably  good  agreement  with  the 
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Re 

cd 

Cpb 

St 

L/d 

Present 

3700 

0.36 

-0.20 

0.22 

90° 

2.60 

104 

0.40 

-0.28 

0.18 

90° 

1.34 

Exp 1 

3700 

0.22 

4200 

-0.23 

104 

0.16 

Exp 2 

3700 

0.21 

104 

0.18 

Comp1  (LES) 

104 

0.393 

0.195 

84°  - 86° 

(DES) 

104 

0.397 

0.200 

O 

L— 

GO 

1 

© 

00 

TABLE  1.  Flow  parameters  of  turbulent  flows  over  a sphere.  Exp1 , 
Kim  & Durbin  (1988);  Exp2,  Sakamoto  & Haniu  (1990);  Comp1, 
Constantinescu  & Squires  (2000).  Here,  Exp  and  Comp  denote  ex- 
perimental and  computational  studies,  respectively. 


previous  results.  Cd  at  Re  = 3700  is  smaller  than  that  at  Re  = 104  in  this 
study,  which  is  closely  related  to  the  fact  that  the  base  pressure  coefficient 
is  larger  at  Re  = 3700.  It  is  interesting  to  note  that  even  though  the  separa- 
tion angle  is  the  same  for  both  Reynolds  numbers,  the  recirculation  bubble 
size  at  Re  = 3700  is  about  twice  that  at  Re  = 104.  This  is  closely  associated 
with  the  characteristics  of  the  separating  shear  layer,  which  is  explained  in 
the  following. 

Figure  2 shows  the  contours  of  the  instantaneous  azimuthal  vorticity  at 
Re  = 3700  and  104  on  the  (a;,  y)-plane.  At  Re  = 3700,  the  separating  shear 
layer  is  elongated  in  the  streamwise  direction  to  form  a cylindrical  vortex 
sheet  and  its  instability  begins  to  appear  at  x « 2d.  On  the  other  hand,  at 
Re  = 104,  the  shear  layer  instability  appears  right  behind  the  sphere.  Note 
that  the  flow  behind  the  sphere  ( x < d)  at  Re  = 3700  contains  few  vortices, 
whereas  at  Re  = 104  many  small-scale  vortices  exist  right  behind  the  sphere. 
Therefore,  at  Re  = 3700,  the  base  pressure  is  higher  and  the  drag  coefficient 
is  smaller  than  at  Re  = 104.  Because  strong  vortical  motion  exists  at  x « 2d 
for  Re  = 3700,  the  lift  fluctuations  are  much  smaller  than  those  at  Re  = 
104.  Due  to  early  onset  of  the  shear  layer  instability  and  strong  vortical 
motions  in  the  near  wake  at  Re  = 104,  the  size  of  the  recirculation  region 
at  Re  = 104  is  much  smaller  than  that  at  Re  = 3700. 

For  clear  understanding  of  the  three-dimensional  structures  in  the  near 
wake,  vortical  structures  are  drawn  in  two  different  ways  in  figure  3,  where 
vortical  surfaces  are  identified  using  the  method  of  Jeong  & Hussain  (1995). 
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(a)  y/d  (c)  y/d 
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(b)  y/d  (d)  y/d 

Figure  4-  Streamwise  velocity  averaged  in  time  and  azimuthal  direction:  (a)  x/d  = 0.2; 

(b)  x/d  = 1.6;  (c)  x/d  = 3.0;  (d)  x/d  = 4.0. , Re  = 3700  (present  study); , 

Re  = 104  (present  study);  •,  Re  = 3700  (Kim  & Durbin  1988). 

That  is,  top  figures  are  the  projection  of  vortical  surfaces  on  the  (x,  y)- 
plane,  whereas  in  bottom  figures  the  vortical  surfaces  are  unfolded  in  the 
azimuthal  direction.  At  Re  = 3700,  the  vortex  sheet  of  the  shear  layer  is 
broken  away  locally  in  the  azimuthal  direction.  However,  at  Re  = 104 , the 
vortex  sheet  is  broken  away  along  the  azimuthal  direction  to  form  vortex 
tubes.  The  term  of  vortex  tube  was  defined  by  Sakamoto  & Haniu  (1990)  as 
the  cylindrical  vortex  shed  from  the  vortex  sheet  separated  from  the  sphere 
surface. 

Figure  4 shows  the  streamwise  velocities  averaged  in  time  and  azimuthal 
direction  at  Re  = 3700  and  104.  The  present  results  at  Re  = 3700  are  in 
good  agreement  with  those  of  single-wire  measurements  by  Kim  & Durbin 
(1988)  at  x/d  = 0.2,  1.6  and  3.0.  A difference  is  observed  at  x = 3.0 d be- 
cause the  measurement  error  of  a single-wire  probe  becomes  large  near  the 
end  of  the  recirculation  region,  where  the  magnitude  of  the  radial  velocity 
is  comparable  to  that  of  the  streamwise  velocity  and  thus  the  flow  angle 
approaching  the  wire  is  not  small.  It  is  clear  that  the  wake  at  Re  = 104 
recovers  more  quickly  than  that  at  Re  = 3700. 

Figures  5 shows  the  contours  of  the  Reynolds  shear  stresses  (— u'xu'r  / u2^) 
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Figure  5.  Contours  of  the  Reynolds  shear  stress  (—u'xu'r/ulo)  in  the  near  wake:  (a)  Re 
= 3700;  (b)  Re  = 104. 


in  the  near  wake  at  Re  = 3700  and  104,  where  the  shear  stress  is  averaged 
in  time  and  azimuthal  direction.  The  Reynolds  shear  stress  as  well  as  the 
normal  stresses  (not  shown  here)  are  zero  up  to  the  separation  point  («s  = 
90°)  for  both  Reynolds  numbers,  showing  that  the  boundary  layer  flow 
before  separation  is  laminar.  It  is  clear  that  the  transition  of  the  wake  to 
turbulence  occurs  right  behind  the  sphere  at  Re  = 104 , whereas  at  Re  = 
3700  the  transition  occurs  far  downstream  (x  ss  2d).  This  phenomenon  is 
closely  related  to  different  onset  of  the  shear  layer  instability  between  two 
Reynolds  numbers. 


4.  Summary 

In  summary,  the  flow  characteristics  between  two  Reynolds  numbers  of 
Re  = 3700  and  104  are  quite  different.  At  Re  = 3700,  the  shear  layer  is 
elongated  in  the  streamwise  direction  to  form  a cylindrical  vortex  sheet 
and  its  instability  begins  to  appear  at  x « 2d.  On  the  other  hand,  at  Re  = 
10  , the  shear  layer  instability  occurs  right  behind  the  sphere  in  the  form 
of  vortex  tube  and  the  flow  becomes  turbulent  in  the  near  wake.  Note  that 
the  flow  behind  the  sphere  (x  < d)  at  Re  = 3700  is  nearly  laminar  and 
contains  few  vortices.  Therefore,  at  Re  = 3700,  the  base  pressure  is  higher 
and  the  drag  coefficient  is  smaller  than  at  Re  = 104.  Because  the  vortices 
exist  at  x ~ 2d  for  Re  = 3700,  the  lift  fluctuations  are  much  smaller  than 
those  at  Re  = 104.  At  Re  = 104,  mixing  in  the  near  wake  is  much  stronger 
than  that  at  Re  = 3700,  and  thus  the  wake  recovers  more  quickly. 
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